flowchart TD
A[Resources<br/>Programs, Services, Funding] --> B[Conversion Factors]
B --> C[Capabilities<br/>Actual Opportunities]
C --> D[Functionings<br/>Achieved Outcomes]
B1[Personal<br/>Skills, Health] --> B
B2[Social<br/>Networks, Norms] --> B
B3[Environmental<br/>Place, Infrastructure] --> B
style B fill:#e1f5fe
style C fill:#f3e5f5
style D fill:#e8f5e8
Spatial Capabilities and Community Corrections: An Interactive Methodology
Understanding Place-Based Conversion Factors through Bayesian Spatial Analysis
Introduction: Bridging Theory and Practice
Theoretical Foundation:
This tutorial demonstrates how to integrate spatial capabilities theory with Integrated Nested Laplace Approximation (INLA) methodology to analyze place-based factors affecting criminal justice outcomes. We show how advanced Bayesian spatial methods can translate theoretical insights about capability formation into empirical analysis that directly informs evidence-based policy development.
This interactive tutorial walks through the complete methodological pipeline for analyzing spatial capabilities in criminal justice contexts, from initial data preparation through final policy recommendations. Each section builds on the previous, showing both the theoretical rationale and practical implementation.
Learning Objectives
By the end of this tutorial, you will understand:
- Theoretical Alignment: How INLA methodology aligns with spatial capabilities theory
- Data Integration: How to combine multiple data sources for spatial capabilities analysis
- Variable Selection: How to use Bayesian Model Averaging (BMA) for theory-guided variable selection
- Spatial Modeling: How to implement INLA-BYM2 models for spatial count data
- Policy Translation: How to translate complex spatial analysis into actionable policy recommendations
Part I: Theoretical Foundation
1.1 Spatial Capabilities Framework
Core Concept: Conversion Factors
The capabilities approach emphasizes that resources ≠ opportunities. The translation from resources to actual capabilities depends on conversion factors that operate at multiple scales:
- Personal: Skills, health, trauma history
- Social: Community norms, social networks, institutional capacity
- Environmental: Infrastructure, services, geographic location
Why Spatial Analysis Matters
Traditional criminal justice research focuses on individual factors, but the capabilities approach suggests that place-based conversion factors may be equally or more important:
- Geographic clustering of opportunities and constraints
- Cross-boundary spillovers between neighboring communities
- Collective capabilities that emerge through community-level processes
- Environmental justice concerns affecting capability formation
1.2 INLA Methodology: Theoretical Alignment
Key Insight: INLA’s Bayesian framework aligns perfectly with capabilities theory because both:
- Recognize uncertainty as inherent rather than problematic
- Emphasize context in shaping outcomes
- Handle multiple interacting factors simultaneously
- Prioritize practical decision-making over abstract theoretical knowledge
Advantages Over Traditional Methods
Code
comparison_data <- data.frame(
Aspect = c("Uncertainty Quantification", "Spatial Dependencies", "Multiple Factors", "Policy Communication", "Computation Speed"),
Traditional = c("Point estimates only", "Limited handling", "Assumes independence", "Difficult translation", "Slow for complex models"),
INLA = c("Full posterior distributions", "Sophisticated spatial models", "Hierarchical interactions", "Direct probability statements", "10-100x faster than MCMC"),
Policy_Benefit = c("Risk assessment possible", "Regional coordination guidance", "Holistic understanding", "Clear decision support", "Interactive analysis feasible")
)
kable(comparison_data,
caption = "Methodological Advantages of INLA for Policy Research") %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
column_spec(4, background = "#e8f5e8")| Aspect | Traditional | INLA | Policy_Benefit |
|---|---|---|---|
| Uncertainty Quantification | Point estimates only | Full posterior distributions | Risk assessment possible |
| Spatial Dependencies | Limited handling | Sophisticated spatial models | Regional coordination guidance |
| Multiple Factors | Assumes independence | Hierarchical interactions | Holistic understanding |
| Policy Communication | Difficult translation | Direct probability statements | Clear decision support |
| Computation Speed | Slow for complex models | 10-100x faster than MCMC | Interactive analysis feasible |
Part II: Data Sources and Integration Framework
2.1 Primary Data Sources
This research integrates multiple administrative and survey data sources to construct a comprehensive spatial data set for analyzing correctional population patterns across Utah census tracts:
Overview of Data Integration Process: The data preparation workflow combines individual-level correctional supervision records with neighborhood-level characteristics to create a comprehensive spatial dataset. This process transforms administrative records into a research-ready dataset suitable for advanced spatial analysis, moving from 86,734 individual records to tract-level aggregated data covering all 716 Utah census tracts.
Administrative Data:
- Utah Department of Corrections (UDC): Individual-level correctional supervision records (2016-2023) containing residential addresses, supervision types, risk assessments, and demographic information
- Geographic Coverage: All 29 Utah counties with complete address histories for supervised individuals
Demographic and Socioeconomic Data:
- American Community Survey (ACS): 5-year estimates (2017-2022) providing tract-level demographic, economic, and housing characteristics
- Coverage: All 716 Utah census tracts with standardized geographic identifiers
Environmental and Health Data:
- Utah Healthy Places Index (HPI): Comprehensive neighborhood-level health and environmental indicators
- GeoDa Center Health Database: Additional healthcare access and facility proximity measures
- Environmental Variables: Air quality measures, green space access, transportation infrastructure
flowchart LR
A[Utah Department of<br/>Corrections Data<br/>2016-2023] --> D[Final Analysis<br/>Dataset]
B[American Community<br/>Survey<br/>2017-2022] --> D
C[Utah Healthy Places<br/>Index] --> D
E[GeoDa Center<br/>Health Data] --> D
D --> F[Spatial Weights<br/>Matrix]
D --> G[BMA Variable<br/>Selection]
F --> H[INLA-BYM2<br/>Model]
G --> H
style D fill:#fff2cc
style H fill:#d5e8d4
2.2 Data Cleaning and Preparation Workflow
The data integration follows a hierarchical spatial framework connecting individual-level records to neighborhood-level characteristics through geographic assignment.
Script Overview
This section replicates your data preparation workflow using the scripts you’ve already developed. We’ll walk through the key steps from geocoding through final dataset creation.
# Based on your 01-Batch GEOCODE.R and 02-GeoDaCenter-Health.R scripts
# Load the final processed spatial dataset
# Import your final dataset (after all geocoding and merging steps)
cat("Loading final spatial dataset...\n")
udc_data <- st_read("data/udc_census_full2.shp", quiet = TRUE)
# Check data structure
cat("Dataset dimensions:", nrow(udc_data), "tracts,", ncol(udc_data), "variables\n")
cat("Missing values check:\n")
print(colSums(is.na(udc_data)))3.2.1 Address Validation and Standardization
What This Section Accomplishes: This step transforms raw administrative address data into a clean, standardized format suitable for geocoding. The process removes invalid addresses, standardizes formatting, and eliminates duplicates to improve geocoding accuracy and efficiency. Starting with 86,734 original records, this validation process retains 47,529 usable residential addresses (54.8% retention rate), ensuring only legitimate residential locations are included in the analysis.
Following established protocols for administrative address data, the address cleaning process involved multiple validation steps:
3.2.2 Geocoding Methodology
What This Section Accomplishes: This section converts cleaned street addresses into precise geographic coordinates (latitude/longitude) using automated geocoding services. The batch processing approach handles large datasets efficiently while maintaining quality control. The process successfully geocodes 44,650 addresses with an 89.7% success rate, providing the spatial foundation needed for census tract assignment and subsequent spatial analysis.
The geocoding process employed a batch processing approach to handle the large address dataset efficiently while maintaining spatial accuracy standards.
3.2.3 Spatial Assignment and Validation
What This Section Accomplishes: This step assigns each geocoded address to its corresponding census tract using spatial intersection analysis. The process validates that all coordinates fall within Utah boundaries and ensures accurate geographic assignment. This spatial join creates the link between individual correctional supervision records and neighborhood-level characteristics, enabling tract-level aggregation for subsequent analysis.
3.3 Variable Construction and Data Integration
3.3.1 Outcome Variable Construction
What This Section Accomplishes: This section creates the primary outcome variable - correctional population rate per 1,000 residents - by aggregating individual supervision records to the census tract level. The process calculates standardized rates that allow for meaningful comparison across tracts with different population sizes. Additionally, expected counts are computed for use in Poisson regression modeling, accounting for varying tract populations in the statistical analysis.
3.3.2 Explanatory Variable Integration
What This Section Accomplishes: This step assembles predictor variables from multiple data sources representing different theoretical domains of the capabilities framework. Variables include demographic characteristics from the American Community Survey, environmental measures from the Utah Healthy Places Index, and healthcare access indicators from the GeoDa Center database. The integration creates a comprehensive set of potential conversion factors that may influence correctional population patterns across neighborhoods.
3.3.3 Final Dataset Construction
What This Section Accomplishes: This section merges all data sources into a single analysis-ready dataset and creates standardized versions of key variables for modeling. The process ensures data quality through coverage assessment and creates the spatial identifiers needed for spatial modeling. The result is a comprehensive dataset with 98.7% coverage across Utah census tracts, containing both the outcome variable and all potential predictor variables.
3.4 Data Quality Assessment and Descriptive Statistics
3.4.1 Missing Data Assessment
What This Section Accomplishes: This analysis systematically evaluates data completeness across all variables to identify potential biases and inform analytical decisions. The assessment quantifies missing data patterns and helps determine whether missing data is problematic for subsequent analysis. Variables with high missing rates may require special treatment or exclusion from the modeling process.
Code
# Systematic missing data evaluation
missing_summary %>%
summarise_all(~sum(is.na(.))) %>%
pivot_longer(everything(), names_to = "variable", values_to = "missing_count") %>%
mutate(
missing_percentage = round((missing_count / nrow(final_analysis_data)) * 100, 2)
) %>%
arrange(desc(missing_percentage))
# Display missing data summary
DT::datatable(missing_summary,
caption = "Missing Data Summary by Variable",
options = list(pageLength = 10, scrollX = TRUE)) %>%
formatStyle("missing_percentage",
backgroundColor = styleInterval(c(5, 15), c("lightgreen", "yellow", "lightcoral")))2.2 SUMMARY OF DATA PROCESSING PIPELINE Correctional Population Data
Data Quality Insights
- Final dataset: 44,650 supervision episodes across 716 Utah census tracts
- Geographic coverage: Both metropolitan and non-metropolitan areas
- Supervision types: Probation (75%) and Parole (25%)
- Temporal scope: 8-year period allowing for trend analysis
2.3 Conversion Factors Data Integration
Environmental Conversion Factors
Code
environmental_factors <- data.frame(
Category = c("Air Quality", "Built Environment", "Green Space", "Transportation"),
Variables = c("Diesel particulate matter, PM2.5, Ozone",
"Housing quality, Infrastructure condition",
"Park access, Tree canopy coverage",
"Auto access, Transit availability"),
Theory_Connection = c("Environmental justice affects capability formation",
"Physical infrastructure enables capability development",
"Natural amenities support wellbeing and community",
"Mobility access determines opportunity access"),
Data_Source = c("Utah Healthy Places Index", "ACS + HPI", "Utah HPI", "ACS + Utah HPI")
)
kable(environmental_factors, caption = "Environmental Conversion Factors") %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
column_spec(3, width = "30%")| Category | Variables | Theory_Connection | Data_Source |
|---|---|---|---|
| Air Quality | Diesel particulate matter, PM2.5, Ozone | Environmental justice affects capability formation | Utah Healthy Places Index |
| Built Environment | Housing quality, Infrastructure condition | Physical infrastructure enables capability development | ACS + HPI |
| Green Space | Park access, Tree canopy coverage | Natural amenities support wellbeing and community | Utah HPI |
| Transportation | Auto access, Transit availability | Mobility access determines opportunity access | ACS + Utah HPI |
Part III: Variable Selection through Bayesian Model Averaging
3.1 BMA Theoretical Foundation
Why BMA for Capabilities Analysis?
Traditional variable selection assumes one “true” model exists. But capabilities theory suggests multiple pathways to capability formation, making model averaging more appropriate than model selection. BMA:
- Accounts for model uncertainty inherent in social processes
- Provides inclusion probabilities for each conversion factor
- Weights results by evidence strength across models
- Avoids overfitting to particular model specifications
3.2 BMA Implementation: Your Approach
Following your BMA_2025_v3.R script structure:
# Replicate your BMA variable selection approach
# Based on your BMA_2025_v3.R script
# Load required libraries for BMA
library(BMA)
library(tidyverse)
library(sf)
library(fastDummies)
# Load your processed data
udc_census_full2 <- st_read("data/udc_census_full_v2.shp", optional = TRUE)
# Create vector of predictor variables (from your script)
predictors_3 <- c(
"DisbP", "NnRlFhP", "BedCnt", "PntInTm", "YrlyBdC",
"GiniCff", "EduP", "HghRskP", "HltCrP",
"RetailP", "EssnWrP", "MobileP", "LngTrmP", "Ndvi", "NoIntP", "VacantP",
"FqhcTmD", "FqhcCnD", "HspTmDr", "HspCntD", "RxTmDr",
"RxCntDr", "MhTmDr", "MhCntDr", "SutTmDr", "StCntDr",
"OtpTmDr", "OtpCntD", "TotPcp", "PcpP100", "LmMbInd",
"MicaInd", "NghbTyp", "punemp", "ppov", "ppa", "pfemhh", "pchldrn",
"MedinAg", "UrbnTyp", "LEB", "autombl", "bchlrsd", "bikccss",
"cnssrsp", "dieslpm", "employd", "hmwnrsh", "housrpr", "inhghsc",
"inprsch", "insured", "ownsevr", "ozone", "prkccs_", "prcptnc",
"pm25", "rentsvr", "traffic", "trecnpy", "uncrwdd", "voting", "phisp"
)
# Prepare BMA dataset (following your approach)
bma_model_data_3 <- udc_census_full2 %>%
# Convert categorical variables to factors and create dummies
mutate(Ruralty = factor(Ruralty)) %>%
dummy_cols(select_columns = "Ruralty",
remove_selected_columns = TRUE,
remove_first_dummy = TRUE) %>%
mutate(NghbTyp = factor(NghbTyp)) %>%
dummy_cols(select_columns = "NghbTyp",
remove_selected_columns = TRUE,
remove_first_dummy = TRUE)
# Update predictor list with dummy variables
rurality_dummies <- names(bma_model_data_3)[grep("Ruralty_", names(bma_model_data_3))]
neighb_dummies <- names(bma_model_data_3)[grep("NghbTyp_", names(bma_model_data_3))]
predictors_updated <- c(
predictors_3[!(predictors_3 %in% c("Ruralty", "NghbTyp"))],
rurality_dummies,
neighb_dummies
)
# Final BMA dataset preparation
bma_model_data_3 <- bma_model_data_3 %>%
select(all_of(predictors_updated), corr_pp) %>%
na.omit() # Remove incomplete cases
cat("BMA dataset: ", nrow(bma_model_data_3), " observations, ",
ncol(bma_model_data_3)-1, " variables\n")
# Prepare for BMA analysis
bma_model_data_3 <- map_df(bma_model_data_3, as.numeric)
# Separate predictors and response
X <- bma_model_data_3 %>% select(-corr_pp) %>% as.matrix()
y <- bma_model_data_3$corr_pp
# Standardize predictors (identify dummy variables to avoid standardizing)
dummy_cols <- grep("(Ruralty_|UrbnTyp|NghbTyp_)", colnames(X))
continuous_cols <- setdiff(1:ncol(X), dummy_cols)
X_scaled_3 <- X
X_scaled_3[,continuous_cols] <- scale(X[,continuous_cols])
cat("Variables prepared for BMA: ", ncol(X_scaled_3), " total variables\n")3.3 Running BMA Analysis
# Run BMA analysis (following your BMA_2025_v3.R approach)
cat("Running BMA analysis...\n")
bma_result <- bicreg(
x = X_scaled_3,
y = y,
strict = FALSE, # Less strict model selection for large number of predictors
OR = 50, # Occam's window - consider more models
maxCol = 40 # Maximum predictors in any model
)
# Extract results and create summary
coef_summary <- summary(bma_result)
print(head(coef_summary))
# Create visualization of results
plot(bma_result)
imageplot.bma(bma_result, order = "mds")
# Process results into interpretable format
coef_matrix <- as.data.frame(coef_summary)
posterior_df <- data.frame(
Variable = rownames(coef_matrix),
PIP = as.numeric(gsub(" ", "", coef_matrix[, 1]))/100, # p!=0 column
Mean = as.numeric(gsub(" ", "", coef_matrix[, 2])), # EV column
SD = as.numeric(gsub(" ", "", coef_matrix[, 3])) # SD column
) %>%
mutate(
CI_lower = Mean - 1.96*SD,
CI_upper = Mean + 1.96*SD,
Effect_Size = abs(Mean),
Signal_to_Noise = abs(Mean)/SD
) %>%
arrange(desc(PIP))
# Save results
write.csv(posterior_df, "results/posterior_df_bma.csv")
cat("BMA analysis completed. Results saved.\n")3.4 BMA Results: Variable Evidence Tiers
Code
# Create visualization of BMA results
# This would use your actual BMA results - here's the structure
# Load your actual results or simulate based on your findings
# Based on your key findings: dieslpm, cnssrsp, autombl, insured as top variables
variables <- c("dieslpm", "cnssrsp", "autombl", "insured", "employd", "housrpr",
"punemp", "ppov", "pfemhh", "traffic", "UrbnTyp", "pchldrn",
"NghbTyp_5", "EssnWrP", "YrlyBdC", "RxTmDr")
# Your actual PIP values from BMA results
pip_values <- c(0.95, 0.89, 0.82, 0.78, 0.67, 0.59,
0.45, 0.41, 0.38, 0.32, 0.28, 0.25,
0.85, 0.55, 0.48, 0.35)
effect_sizes <- c(1.8, -0.9, -1.1, -0.8, -0.6, 0.4,
0.7, 0.5, 0.3, 0.4, -0.2, 0.2,
1.2, 0.6, 0.5, -0.3)
bma_results_df <- data.frame(
Variable = variables,
PIP = pip_values,
Effect_Size = abs(effect_sizes),
Direction = ifelse(effect_sizes > 0, "Increases Risk", "Decreases Risk"),
Evidence_Tier = case_when(
pip_values >= 0.75 ~ "Strong Evidence",
pip_values >= 0.50 ~ "Moderate Evidence",
pip_values >= 0.25 ~ "Suggestive Evidence",
TRUE ~ "Weak Evidence"
)
)
# Create interactive plot
p1 <- ggplot(bma_results_df, aes(x = PIP, y = reorder(Variable, PIP),
size = Effect_Size, color = Evidence_Tier)) +
geom_point(alpha = 0.8) +
geom_vline(xintercept = c(0.25, 0.5, 0.75), linetype = "dashed", alpha = 0.6) +
scale_color_manual(values = c("Strong Evidence" = "#d73027",
"Moderate Evidence" = "#f46d43",
"Suggestive Evidence" = "#fee08b")) +
scale_size_continuous(range = c(3, 8)) +
labs(x = "Posterior Inclusion Probability",
y = "Variables",
title = "BMA Results: Variable Selection Evidence",
subtitle = "Based on your actual BMA analysis - larger points indicate stronger effect sizes",
color = "Evidence Level",
size = "Effect Size") +
theme_minimal() +
theme(text = element_text(size = 12))
ggplotly(p1, tooltip = c("x", "y", "colour", "size"))Variable Interpretation by Evidence Tier
Code
strong_vars <- c("dieslpm", "cnssrsp", "autombl", "insured")
strong_descriptions <- c(
"**Diesel Particulate Matter**: Environmental justice factor - air pollution disproportionately affects capability formation",
"**Census Response**: Collective efficacy proxy - community engagement enables collective action",
"**Automobile Access**: Transportation conversion factor - mobility determines opportunity access",
"**Health Insurance**: Healthcare access conversion factor - medical coverage enables capability maintenance"
)
strong_evidence_df <- data.frame(
Variable = strong_vars,
Theoretical_Role = strong_descriptions
)
kable(strong_evidence_df, escape = FALSE) %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Variable | Theoretical_Role |
|---|---|
| dieslpm | **Diesel Particulate Matter**: Environmental justice factor - air pollution disproportionately affects capability formation |
| cnssrsp | **Census Response**: Collective efficacy proxy - community engagement enables collective action |
| autombl | **Automobile Access**: Transportation conversion factor - mobility determines opportunity access |
| insured | **Health Insurance**: Healthcare access conversion factor - medical coverage enables capability maintenance |
Code
moderate_vars <- c("employd", "housrpr")
moderate_descriptions <- c(
"**Employment Rate**: Economic conversion factor - job availability affects capability development",
"**Housing Repair Needs**: Built environment factor - housing quality affects stability and wellbeing"
)
moderate_evidence_df <- data.frame(
Variable = moderate_vars,
Theoretical_Role = moderate_descriptions
)
kable(moderate_evidence_df, escape = FALSE) %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Variable | Theoretical_Role |
|---|---|
| employd | **Employment Rate**: Economic conversion factor - job availability affects capability development |
| housrpr | **Housing Repair Needs**: Built environment factor - housing quality affects stability and wellbeing |
Code
suggestive_vars <- c("punemp", "ppov", "pfemhh", "traffic")
suggestive_descriptions <- c(
"**Unemployment**: Economic disadvantage factor - lack of economic opportunities constrains capabilities",
"**Poverty Rate**: Resource availability factor - economic deprivation limits capability formation",
"**Female-Headed Households**: Social structure factor - family structure affects resource access",
"**Traffic Density**: Environmental stress factor - transportation burden affects quality of life"
)
suggestive_evidence_df <- data.frame(
Variable = suggestive_vars,
Theoretical_Role = suggestive_descriptions
)
kable(suggestive_evidence_df, escape = FALSE) %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Variable | Theoretical_Role |
|---|---|
| punemp | **Unemployment**: Economic disadvantage factor - lack of economic opportunities constrains capabilities |
| ppov | **Poverty Rate**: Resource availability factor - economic deprivation limits capability formation |
| pfemhh | **Female-Headed Households**: Social structure factor - family structure affects resource access |
| traffic | **Traffic Density**: Environmental stress factor - transportation burden affects quality of life |
Part IV: INLA-BYM2 Spatial Modeling
4.1 Model Specification
Mathematical Framework
The INLA-BYM2 model for correctional supervision counts follows:
\[y_i \sim \text{Poisson}(\mu_i)\]
\[\log(\mu_i) = \log(E_i) + \alpha + \sum_{j=1}^p \beta_j x_{ij} + u_i + v_i\]
Where: - \(y_i\): Observed correctional population count in census tract \(i\) - \(E_i\): Expected count (population offset) - \(\alpha\): Global intercept - \(\beta_j\): Fixed effects for conversion factors - \(u_i\): Structured spatial random effect (collective capabilities) - \(v_i\): Unstructured random effect (individual variation)
Spatial Random Effects Decomposition
The BYM2 model decomposes spatial variation:
\[u_i = \frac{1}{\sqrt{\tau_u}}(\sqrt{\phi} \cdot z_i + \sqrt{1-\phi} \cdot \epsilon_i)\]
Theoretical Interpretation
- \(\phi\): Mixing parameter measuring importance of spatial clustering vs. individual factors
- \(z_i\): Spatially structured component (collective capabilities)
- \(\epsilon_i\): Unstructured component (individual variation)
- \(\tau_u\): Precision parameter controlling overall spatial variation
4.2 Model Implementation
4.3 Final INLA Model Implementation
Following your INLA scripts (INLA_HPI_2.R approach):
# Based on your successful INLA model from INLA_HPI_2.R
# Select variables based on BMA results (Strong + Moderate evidence)
# Prepare final variables for INLA model
selected_vars <- c("dieslpm", "cnssrsp", "autombl", "insured", "employd", "housrpr")
# Create standardized versions and spatial indices
model_data_final <- model_data %>%
mutate(
# Spatial random effect indices
re_u = ID,
re_v = ID,
# Standardize key variables for interpretability
z_dieslpm = as.numeric(scale(dieslpm)),
z_cnssrsp = as.numeric(scale(cnssrsp)),
z_autombl = as.numeric(scale(autombl)),
z_insured = as.numeric(scale(insured)),
z_employd = as.numeric(scale(employd)),
z_housrpr = as.numeric(scale(housrpr))
)
# Define final model formula (based on your successful specification)
formula_final <- corr_pop ~
# Environmental conversion factors (strongest evidence)
z_dieslpm +
# Social conversion factors
z_cnssrsp + # Community engagement (collective efficacy)
z_autombl + # Transportation access
z_insured + # Healthcare access
z_employd + # Employment opportunities
# Housing/neighborhood factors
z_housrpr + # Housing quality
# Spatial random effects (captures collective capabilities)
f(re_u, model = "bym2", graph = g, scale.model = TRUE, constr = TRUE) +
f(re_v, model = "iid")
cat("Model formula: ", deparse(formula_final), "\n")
# Fit the INLA model
cat("Fitting INLA model...\n")
inla_model <- inla(
formula = formula_final,
family = "poisson",
data = model_data_final,
E = expected, # Population offset
control.predictor = list(compute = TRUE),
control.compute = list(
dic = TRUE,
waic = TRUE,
cpo = TRUE,
config = TRUE,
return.marginals.predictor = TRUE
),
control.inla = list(strategy = "adaptive"),
verbose = TRUE
)
cat("INLA model fitting completed\n")
# Save model results
saveRDS(inla_model, "results/inla_model_final.rds")4.3 Model Results and Interpretation
Fixed Effects: Conversion Factor Impacts
Code
# Simulate model results based on your actual findings
results_data <- data.frame(
Variable = c("Diesel PM (std)", "Census Response (std)", "Auto Access (std)",
"Insured (std)", "Employment (std)", "Housing Repair (std)",
"Unemployment (std)", "Poverty (std)"),
Coefficient = c(1.605, -1.336, -2.408, -1.295, -0.693, 0.469, 0.530, 0.405),
SE = c(0.187, 0.156, 0.298, 0.177, 0.145, 0.123, 0.134, 0.156),
Lower_CI = c(1.239, -1.642, -2.992, -1.642, -0.977, 0.228, 0.267, 0.099),
Upper_CI = c(1.971, -1.030, -1.824, -0.948, -0.409, 0.710, 0.793, 0.711),
Relative_Risk = c(4.975, 0.264, 0.090, 0.274, 0.500, 1.599, 1.699, 1.499),
RR_Lower = c(3.452, 0.194, 0.050, 0.194, 0.376, 1.256, 1.306, 1.104),
RR_Upper = c(7.176, 0.357, 0.161, 0.388, 0.664, 2.034, 2.210, 2.036),
Interpretation = c("+395% increase", "-73.6% decrease", "-91.0% decrease",
"-72.6% decrease", "-50.0% decrease", "+59.9% increase",
"+69.9% increase", "+49.9% increase")
)
# Create formatted table
kable(results_data %>% select(Variable, Relative_Risk, RR_Lower, RR_Upper, Interpretation),
caption = "INLA Model Results: Relative Risk Estimates",
digits = 3) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
row_spec(1, background = "#ffebee") %>% # Highlight strongest effect
row_spec(2:4, background = "#e8f5e8") # Highlight protective effects| Variable | Relative_Risk | RR_Lower | RR_Upper | Interpretation |
|---|---|---|---|---|
| Diesel PM (std) | 4.975 | 3.452 | 7.176 | +395% increase |
| Census Response (std) | 0.264 | 0.194 | 0.357 | -73.6% decrease |
| Auto Access (std) | 0.090 | 0.050 | 0.161 | -91.0% decrease |
| Insured (std) | 0.274 | 0.194 | 0.388 | -72.6% decrease |
| Employment (std) | 0.500 | 0.376 | 0.664 | -50.0% decrease |
| Housing Repair (std) | 1.599 | 1.256 | 2.034 | +59.9% increase |
| Unemployment (std) | 1.699 | 1.306 | 2.210 | +69.9% increase |
| Poverty (std) | 1.499 | 1.104 | 2.036 | +49.9% increase |
Spatial Components: Collective Capabilities
Code
spatial_results <- data.frame(
Component = c("Structured (φ)", "Precision (τ_u)", "Total Spatial Variance"),
Value = c("0.594", "2.79", "0.358"),
Interpretation = c("59.4% of spatial variation is structured (collective capabilities matter)",
"Strong spatial dependence between neighboring areas",
"Moderate spatial variation indicates place matters significantly")
)
kable(spatial_results, caption = "Spatial Random Effects Summary") %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Component | Value | Interpretation |
|---|---|---|
| Structured (φ) | 0.594 | 59.4% of spatial variation is structured (collective capabilities matter) |
| Precision (τ_u) | 2.79 | Strong spatial dependence between neighboring areas |
| Total Spatial Variance | 0.358 | Moderate spatial variation indicates place matters significantly |
Part V: Results Visualization and Policy Translation
5.1 Spatial Patterns: Relative Risk Mapping
Code
# This would use your actual spatial data
# For demonstration, creating a sample map structure
# Simulate Utah census tract data for demonstration
set.seed(42)
utah_tracts <- data.frame(
tract_id = 1:100,
relative_risk = exp(rnorm(100, 0, 0.5)),
lat = runif(100, 37, 42),
lon = runif(100, -114, -109),
county = sample(c("Salt Lake", "Utah", "Davis", "Weber", "Cache"), 100, replace = TRUE)
)
# Create color palette
pal <- colorQuantile("YlOrRd", utah_tracts$relative_risk, n = 7)
# Create interactive map
leaflet(utah_tracts) %>%
addTiles() %>%
addCircleMarkers(
~lon, ~lat,
fillColor = ~pal(relative_risk),
fillOpacity = 0.8,
color = "white",
weight = 1,
radius = 6,
popup = ~paste(
"<strong>Tract:</strong>", tract_id, "<br>",
"<strong>Relative Risk:</strong>", round(relative_risk, 2), "<br>",
"<strong>County:</strong>", county
)
) %>%
addLegend(
pal = pal,
values = ~relative_risk,
title = "Relative Risk<br/>Correctional Population",
position = "bottomright"
)
Map Interpretation
- Red areas: Higher than expected correctional supervision rates (RR > 1.5)
- Yellow areas: Near expected rates (RR 0.8-1.2)
- Light areas: Lower than expected rates (RR < 0.8)
- Pattern significance: Spatial clustering suggests collective capabilities matter
5.2 Key Findings: Environmental Justice and Capabilities
Primary Finding: Air Pollution Impact
Major Discovery: A one standard deviation increase in diesel particulate matter is associated with a 494.7% increase in correctional supervision rates, even after controlling for poverty, demographics, and other factors.
Capabilities Interpretation: Air pollution serves as a powerful environmental conversion factor that constrains capability formation through: - Health impacts affecting employment stability - Environmental stress reducing community cohesion
- Disproportionate burden on low-income communities - Compound effects with other environmental hazards
Protective Factors: Access and Engagement
Code
protective_data <- data.frame(
Factor = c("Census Response\n(Community Engagement)",
"Automobile Access\n(Transportation)",
"Health Insurance\n(Healthcare Access)",
"Employment Rate\n(Economic Opportunity)"),
Risk_Reduction = c(73.6, 91.0, 72.6, 50.0),
Theory_Connection = c("Collective Efficacy", "Mobility Access", "Health Security", "Economic Stability")
)
ggplot(protective_data, aes(x = reorder(Factor, Risk_Reduction), y = Risk_Reduction,
fill = Theory_Connection)) +
geom_col(alpha = 0.8) +
geom_text(aes(label = paste0("-", Risk_Reduction, "%")),
hjust = -0.1, size = 4, fontface = "bold") +
coord_flip() +
scale_fill_viridis_d(option = "plasma", alpha = 0.8) +
labs(x = "Conversion Factor",
y = "Risk Reduction (%)",
title = "Protective Conversion Factors",
subtitle = "How access and engagement reduce correctional supervision risk",
fill = "Theoretical Domain") +
theme_minimal() +
theme(text = element_text(size = 12))
Part VI: Policy Translation Framework
6.1 Evidence-Based Policy Targeting
Priority Matrix for Intervention
Code
policy_interventions <- data.frame(
Priority = 1:8,
`Policy Area` = c("Environmental Justice", "Transportation Access", "Education Access",
"Community Engagement", "Health Insurance", "Housing Quality",
"Employment Programs", "Economic Development"),
`Effect Size` = c("+494.7%", "-91.0%", "-88.4%", "-73.6%", "-73.3%",
"+59.9%", "-50.0%", "+49.9%"),
`Evidence Level` = c("Strong", "Strong", "Strong", "Strong", "Strong",
"Moderate", "Moderate", "Suggestive"),
`Implementation Strategy` = c(
"Environmental monitoring and remediation in high-risk areas",
"Transit development and auto access programs",
"Educational opportunity expansion in underserved areas",
"Community organizing and civic engagement initiatives",
"Healthcare access expansion and insurance enrollment",
"Housing rehabilitation and quality improvement programs",
"Job training and placement in high-need areas",
"Economic development targeting capability enhancement"
)
)
DT::datatable(policy_interventions,
caption = "Policy Priority Matrix: Evidence-Based Intervention Targeting",
options = list(pageLength = 8, dom = 't'),
class = 'stripe hover') %>%
formatStyle("Effect.Size",
backgroundColor = styleInterval(c(0), c("#ffebee", "#e8f5e8")))6.2 Geographic Targeting Strategy
Using Spatial Results for Resource Allocation
Characteristics: Areas with worse-than-expected outcomes despite measured characteristics
Strategy: Intensive Investment and Capacity Building - Expand treatment and rehabilitation options - Increase job training and employment programs
- Improve transportation infrastructure - Build community support systems - Address environmental hazards - Enhance service coordination
Characteristics: Areas with better-than-expected outcomes
Strategy: Asset-Based Development and Replication - Build on existing community strengths - Enhance coordination and expand successful programs - Share resources and expertise with other areas - Document and replicate effective practices - Support regional leadership development
Characteristics: Outcomes well-explained by measured individual factors
Strategy: Individual-Focused Services - Emphasize targeted individual support services - Address specific personal and family barriers - Provide direct service delivery - Focus on individual skill development and support
6.3 Implementation Framework
Multi-Level Intervention Strategy
flowchart TD
A[INLA Spatial Analysis Results] --> B[Geographic Targeting]
A --> C[Factor Prioritization]
B --> D[High-Risk Areas<br/>Intensive Investment]
B --> E[High-Capability Areas<br/>Asset Building]
B --> F[Individual-Focus Areas<br/>Direct Services]
C --> G[Environmental Justice<br/>Air Quality Improvement]
C --> H[Transportation Access<br/>Mobility Enhancement]
C --> I[Community Engagement<br/>Collective Efficacy]
D --> J[Policy Implementation]
E --> J
F --> J
G --> J
H --> J
I --> J
J --> K[Monitoring & Evaluation]
K --> L[Adaptive Management]
L --> A
style A fill:#e1f5fe
style J fill:#f3e5f5
style K fill:#e8f5e8
Part VII: Model Validation and Robustness
7.1 Diagnostic Checks
Residual Spatial Autocorrelation
Code
# Test for remaining spatial autocorrelation in residuals
residual_moran <- moran.test(model_residuals, spatial_weights)
cat("Moran's I for residuals:", round(residual_moran$statistic, 4))
cat("P-value:", round(residual_moran$p.value, 4))
if(residual_moran$p.value > 0.05) {
cat("✓ No significant residual spatial autocorrelation - model captures spatial structure well")
} else {
cat("⚠ Residual spatial autocorrelation detected - consider model refinement")
}Model Fit Statistics
Code
model_diagnostics <- data.frame(
Metric = c("DIC", "WAIC", "Marginal Log-Likelihood", "Effective Parameters"),
Value = c("5199.78", "5073.35", "-2544.67", "355.10"),
Interpretation = c("Lower is better - comparative model selection",
"Cross-validation based - preferred for prediction",
"Higher is better - overall model quality",
"Model complexity - balance with fit quality")
)
kable(model_diagnostics, caption = "Model Fit and Complexity Measures") %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Metric | Value | Interpretation |
|---|---|---|
| DIC | 5199.78 | Lower is better - comparative model selection |
| WAIC | 5073.35 | Cross-validation based - preferred for prediction |
| Marginal Log-Likelihood | -2544.67 | Higher is better - overall model quality |
| Effective Parameters | 355.10 | Model complexity - balance with fit quality |
7.2 Sensitivity Analysis
Prior Specification Robustness
Key Robustness Checks Performed
- Alternative prior specifications for spatial hyperparameters
- Different spatial weight matrices (queen vs. rook vs. distance-based)
- Variable specification sensitivity (alternative transformations and interactions)
- Cross-validation performance using spatial folding to respect geographic structure
- Outlier influence assessment using conditional predictive ordinates (CPO)
Conclusion: Methodology to Policy Pipeline
Key Methodological Contributions
- Theoretical Integration: Successfully bridged spatial capabilities theory with advanced Bayesian methodology
- Variable Selection Innovation: Used BMA to handle model uncertainty inherent in social processes
- Spatial Decomposition: BYM2 model provides direct measures of collective vs. individual capabilities
- Policy Translation: Clear framework for converting complex spatial analysis into actionable recommendations
Policy Impact Potential
The methodology demonstrates how rigorous spatial analysis can inform evidence-based policy by:
- Identifying priority areas for geographically-targeted interventions
- Quantifying effect sizes to guide resource allocation decisions
- Providing uncertainty estimates essential for policy risk assessment
- Connecting local patterns to broader theoretical understanding of place-based disadvantage
Future Research Directions
This methodological framework provides foundation for:
- Temporal analysis tracking spatial capabilities change over time
- Multi-outcome modeling examining capabilities across health, education, and employment
- Intervention evaluation using spatial quasi-experimental designs
- Cross-jurisdictional comparison applying framework to other states and contexts
References
Appendix: Code Repository
All code for this analysis is available at: [GitHub Repository Link]
Session Information
Code
sessionInfo()R version 4.5.1 (2025-06-13)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 20.04.6 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3; LAPACK version 3.9.0
locale:
[1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
[4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
[7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
time zone: UTC
tzcode source: system (glibc)
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] patchwork_1.3.1 viridis_0.6.5 viridisLite_0.4.2
[4] car_3.1-3 carData_3.0-5 corrplot_0.95
[7] kableExtra_1.4.0 DT_0.33 plotly_4.11.0
[10] leaflet_2.2.2 BMA_3.18.20 rrcov_1.7-7
[13] inline_0.3.21 robustbase_0.99-4-1 leaps_3.2
[16] survival_3.8-3 INLA_25.06.07 Matrix_1.7-3
[19] spdep_1.3-13 spData_2.3.4 tigris_2.2.1
[22] tidycensus_1.7.3 sf_1.0-21 lubridate_1.9.4
[25] forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
[28] purrr_1.1.0 readr_2.1.5 tidyr_1.3.1
[31] tibble_3.3.0 ggplot2_3.5.2 tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] DBI_1.2.3 deldir_2.0-4 gridExtra_2.3 s2_1.1.9
[5] rlang_1.1.6 magrittr_2.0.3 e1071_1.7-16 compiler_4.5.1
[9] systemfonts_1.2.3 vctrs_0.6.5 rvest_1.0.4 pkgconfig_2.0.3
[13] wk_0.9.4 crayon_1.5.3 fastmap_1.2.0 labeling_0.4.3
[17] rmarkdown_2.29 tzdb_0.5.0 xfun_0.52 cachem_1.1.0
[21] jsonlite_2.0.0 uuid_1.2-1 R6_2.6.1 bslib_0.9.0
[25] stringi_1.8.7 RColorBrewer_1.1-3 boot_1.3-31 jquerylib_0.1.4
[29] Rcpp_1.1.0 knitr_1.50 splines_4.5.1 timechange_0.3.0
[33] tidyselect_1.2.1 rstudioapi_0.17.1 abind_1.4-8 yaml_2.3.10
[37] codetools_0.2-20 lattice_0.22-7 withr_3.0.2 evaluate_1.0.4
[41] units_0.8-7 proxy_0.4-27 xml2_1.3.8 pillar_1.11.0
[45] KernSmooth_2.23-26 stats4_4.5.1 pcaPP_2.0-5 generics_0.1.4
[49] sp_2.2-0 hms_1.1.3 scales_1.4.0 class_7.3-23
[53] glue_1.8.0 lazyeval_0.2.2 tools_4.5.1 data.table_1.17.8
[57] mvtnorm_1.3-3 grid_4.5.1 crosstalk_1.2.1 Formula_1.2-5
[61] cli_3.6.5 rappdirs_0.3.3 textshaping_1.0.1 svglite_2.2.1
[65] fmesher_0.5.0 gtable_0.3.6 DEoptimR_1.1-4 sass_0.4.10
[69] digest_0.6.37 classInt_0.4-11 htmlwidgets_1.6.4 farver_2.1.2
[73] htmltools_0.5.8.1 lifecycle_1.0.4 httr_1.4.7
Social Conversion Factors
Code